
#include <fstream>
#include "mesh.h"

/// A functor for computing plane equations
struct PlaneEquation {
   template <class Facet>
   typename Facet::Plane_3 operator()(Facet& f) {
      typename Facet::Halfedge_handle h = f.halfedge();
      typedef typename Facet::Plane_3  Plane;
      return Plane(h->vertex()->point(),
                   h->next()->vertex()->point(),
                   h->next()->next()->vertex()->point());
   }
};

/// A functor for computing vertex normals 
struct VertexNormal {
   template <class Vertex>
   void operator()(Vertex& v) {
      Mesh::Halfedge_around_vertex_const_circulator 
            hbegin = v.vertex_begin(), h = hbegin;
      int n = 0;
      Vector3 sum(0,0,0);
      do {
         if (not h->is_border()) {
            sum = sum + h->facet()->plane().orthogonal_vector();
            n++;
         }
      } while (++h != hbegin);
      if (n>0 and sum*sum>0) {
         sum = sum/n;
         sum = sum/sqrt(sum*sum);
      }
      v.normal(sum);
   }
};

Mesh::Mesh() {
   
   _nvertices=_nfaces=0;
}

Mesh::Mesh(const char* filename) {
   
   _nvertices=_nfaces=0;
   loadShape(filename);
}

Mesh::~Mesh() {
   
}

void Mesh::loadShape(const char *filename) {
   
   cout<<"Loading... ";
   
   std::ifstream fin(filename);
   fin>>*this;
   fin.close();
   auxiliarData();
}

void Mesh::saveShapeOFF(const char *filename) {
   
   if (_nvertices==0 or _nfaces==0)   return;
   
   std::ofstream fout(filename);
   fout<<*this;
   fout.close();
}

Mesh* Mesh::load_off(const char *filename) {
   
   std::fstream off_file(filename, std::ios_base::in);
   if (off_file) {
      Mesh *p = new Mesh(filename);
      return p;
   }
   return 0; 
}

void Mesh::auxiliarData() {
   
   if (not this->is_pure_triangle()) {
      cout<<"### *** Attention . Attention . Attention . Attention *** ###"<<endl;
      cout<<"This example requires the mesh model to have triangle faces only!!!!."<<endl;
      return;
   }
   
   this->normalize_border();
   _nvertices = this->size_of_vertices();
   _nfaces    = this->size_of_facets();
   
   cout<<endl<<"Loading Done!"<<endl;
   cout<<"Vertices Number: "<<_nvertices<<endl;
   cout<<"Triangles Number: "<<_nfaces<<endl;
   if (this->size_of_border_edges() > 0)
      cout<<"The mesh has border edges "<<endl;
   else
      cout<<"The mesh hasnt border edges "<<endl;
   cout<<"KBytes used by the model: "<<this->bytes()/1000<<endl;
   cout<<"KBytes reserved by the model: "<<this->bytes_reserved()/1000<<endl;
   
   /// Compute the normal vector of each face
   std::transform(this->facets_begin(), this->facets_end(), 
                  this->planes_begin(), PlaneEquation());
   
   /// Compute the normal vector of each vertex
   std::for_each(this->vertices_begin(), this->vertices_end(), VertexNormal());
   
   /// Build the array containing the index of faces
   buildFacesIndices();
   
   /// Computes the limits of the model
   limits();
   
   /// Centers and scales the model
   //~ centerAndScale();
}

void Mesh::buildFacesIndices() {
   
   cout<<"buildArrays begin ... "<<endl;
   faces.resize(3*this->size_of_facets());
   // Build the array of faces
   unsigned iif=0;
   for (Facet_iterator fi=this->facets_begin(); fi!=this->facets_end(); ++fi, iif+=3) {
      Halfedge_around_facet_circulator j = fi->facet_begin();
      
      const unsigned &id0 = std::distance(this->vertices_begin(), j->vertex());   ++j;
      const unsigned &id1 = std::distance(this->vertices_begin(), j->vertex());   ++j;
      const unsigned &id2 = std::distance(this->vertices_begin(), j->vertex());
      faces[iif  ]=id0;
      faces[iif+1]=id1;
      faces[iif+2]=id2;
   }
   cout<<"... done! "<<endl;
}

void Mesh::limits() {
   
   Vertex_iterator vit = this->vertices_begin();
   const Point3 &p = vit->point();
   for (int i = 0; i < 3; ++i) { mincoo[i] = maxcoo [i] = p[i]; }
   
   while (vit != this->vertices_end()) {
      const Point3 &p = vit->point();
      for (int i = 0; i < 3; ++i) {
         if (p[i] < mincoo[i])      mincoo[i] = p[i];
         else if (p[i] > maxcoo[i]) maxcoo[i] = p[i];
      }
      ++vit;
   }
   
   /// building the minimum bounding cube 
   Point3 WorldpMin(mincoo[0],mincoo[1],mincoo[2]);
   Point3 WorldpMax(maxcoo[0],maxcoo[1],maxcoo[2]);
   
   cout<<"pmin: "<<mincoo[0]<<" "<<mincoo[1]<<" "<<mincoo[2]<<endl;
   cout<<"pmax: "<<maxcoo[0]<<" "<<maxcoo[1]<<" "<<maxcoo[2]<<endl;
   
   Vector3 size(WorldpMin, WorldpMax);
   size = size*0.5;
   Point3 center = WorldpMin + size;
   double radius = sqrt (size * size); // vector lenght
   size = Vector3(radius, radius, radius);
   WorldpMin = center - size;
   WorldpMax = center + size;
   
   mincoo[0]=WorldpMin[0]; mincoo[1]=WorldpMin[1]; mincoo[2]=WorldpMin[2];
   maxcoo[0]=WorldpMax[0]; maxcoo[1]=WorldpMax[1]; maxcoo[2]=WorldpMax[2];
}

void Mesh::centerAndScale(Number ScaleFactor) {
   
   assert(ScaleFactor > 0);
   
   Number dimv[]={maxcoo[0]-mincoo[0], maxcoo[1]-mincoo[1], maxcoo[2]-mincoo[2]};
   assert(dimv[0]>0 and dimv[1]>0 and dimv[2]>0);
   
   // Determine the scale factor
   Number Scale;
   if (dimv[0] >= dimv[1] && dimv[0] >= dimv[2]) Scale = ScaleFactor/dimv[0];
   else if (dimv[1] >= dimv[0] && dimv[1] >= dimv[2]) Scale = ScaleFactor/dimv[1];
   else Scale = ScaleFactor/dimv[2];
   
   // Determine the translation vector
   Number transv[]={mincoo[0]+maxcoo[0], mincoo[1]+maxcoo[1], mincoo[2]+maxcoo[2]};
   transv[0] *= 0.5f;
   transv[1] *= 0.5f;
   transv[2] *= 0.5f;
   
   Vertex_iterator vit = this->vertices_begin();
   while (vit != this->vertices_end()) {
      Point3 &p = vit->point();
      p = Point3(p[0]-transv[0], p[1]-transv[1], p[2]-transv[2]);
      p = Point3(p[0]*Scale, p[1]*Scale, p[2]*Scale);
      ++vit;
   }
}
